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The t-J model is analysed in the limit of strong anisotropy, where the transverse components 
of electron spin are neglected. We propose a slave—particle-type approach that is valid, in 
contradiction to many of the standard approaches, in the low—doping regime and becomes ex¬ 
act for a half-filled system. We describe an effective method that allows to numerically study 
the system with the no-double—occupancy constraint rigorously taken into account at each 
lattice site. Then, we use this approach to demonstrate the destruction of the antiferromag¬ 
netic order by increasing doping and formation of Nagaoka polarons in the strong interaction 
regime. 
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1. Introduction 

It is commonly believed that the richness of the behaviour of strongly correlated 
systems is a result of a competition between the kinetic and interaction energies 
[T]. Unfortunately, due to the presence of strong correlations many of the ’’tradi¬ 
tional” solid state methods, like the density functional theory within local density 
approximation, or many-body perturbation theory, that handled impressively well 
simple metals, covalent semiconductors, closed-shell ionic insulators, and even in- 
termetallic compounds, cannot be used. It has been recognized for many years that 
strongly-correlated systems require a distinct paradigm from what was successful 
for the mentioned above systems. 

It is also believed that the essence of the physics of the strongly correlated systems 
can be described by simple one-band Hamiltonians that are able to properly take 
into account the competition between the kinetic and interaction energies. Two of 
the most acceptable models are the Hubbard model [2] and its effective strong- 
interaction version, namely the t-J model 13 HI- Both these models contributed 
greatly to our understanding of strongly correlated systems. Unfortunately, apart 
from some specific cases, none of these models can be solved exactly. Therefore, it is 
very important to develop analytical or numerical methods that can be applied to 
systems described by interacting Hamiltonians. Moreover, it is equally important 
to be able to determine the errors introduced by the applied approximation. 
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1.1. The Hubbard and t-J models 

The Hamiltonian of the Hubbard model, originally introduced to describe correla¬ 
tion effects in narrow d-band materials, has the following form: 

-f^Hubb — t ^ ^ ^ ^ ^ (1) 

{ij)a i 


where the first term describes the kinetic energy and the second the interactions. 
Here, c\^ (cia) creates (annihilates) an electron of spin a at site i and the occupation 
number operator Uia = Hilbert space of the Hubbard model contains 

four states per site: |0), | t); | 1) and | tl)- 

Since in the Hubbard model there is only an on-site interaction, in the limit of 
large U it is energetically very expensive for electrons to hop onto already occupied 
sites. Therefore, for the average occupation less or equal to one electron per lattice 
site the low energy processes take place mainly in the lower Hubbard subband. 
However, virtual excitations with double occupied sites may increase the electron 
mobility leading to lowering the total energy. The effective Hamiltonian can be 
derived from the strong coupling expansion of the Hubbard model with respect 
to t/U. It was shown that that virtual excitations generates a spin-spin exchange 
interaction between neighbouring sites, the so-called kinetic exchange. After the 
transformation the states in the lower Hubbard band are described by the t-J 
model acting in a projected Hilbert space containing only three states per site: 
|0), I t), I 4 ,). The state | ti) is removed by the Gutzwiller projection operator. 
The Hamiltonian of the t-J model is given by: 

Ht-j = —t (^i^j ~ > (2) 

(bV (b> 

where the antiferromagnetic exchange constant J = At^/U. Cia (c|^) represents 
fermionic annihilation (creation) operators projected onto a space without dou¬ 
ble occupancy: = (1 — ni-o-)cio-- Despite a potential inadequacy of this model 

to represent real strongly correlated materials, it is still the simplest model that 
captures the important antiferromagnetic correlations of weakly doped antiferro- 
magnets. Thus, it is crucial that the properties of this model are well understood. 
The Hamiltonian given by Eq. Q has been investigated intensively by different an¬ 
alytical and numerical methods. The analytical methods are usually limited to only 
one or two holes in an antiferromagnetic background. It is very difficult to treat in 
a systematic non-perturbative way systems with strong correlations. In the case 
of the t-J model an additional difficulty comes from the fact that the operators 
Cirr and c|^ do not fulfil the usual fermionic commutation rules. This non-fermionic 
behaviour results, in turn, from the projection of the states with doubly occu¬ 
pied lattice sites. Unfortunately, the constraint of no double occupancy becomes 
very important close to half filling and only methods which are capable of taking it 
into account without uncontrollable approximations can give reliable results in this 
regime. And this is a regime of particular interest because the high-temperature 
superconductors are slightly doped antiferromagnets. 

Due to the difficulties in analytical approaches, numerical methods such as exact 
diagonalization, density-matrix renormalization group and quantum Monte-Carlo 
are extensively performed to study this model. The exact diagonalization can only 
be performed in a very small lattice size, and the density-matrix renormalization 
group method is largely restricted to one-dimensional systems. In contrast, Quan- 
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turn Monte Carlo simulation is the only systematic and scalable method with suf¬ 
ficient numerical accuracy for higher dimensional problems. However, this method 
also has the notorious fermion sign problem which makes low temperature proper¬ 
ties inaccessible. 


1.2. Slave particle approaches to the t-J model 

The single occupancy constraint, that makes analytical approaches to the t-J 
model so difficult can be written as 

X!/<1, (3) 


for every lattice site i. In order to treat this constraint in a controllable way a 
number of slave-particle methods have been proposed [SHa]. In the slave particle 
formalism, the electron operator is expressed in terms of auxiliary fermions and 
bosons. For instance, in the slave boson formalism the electron annihilation oper¬ 
ator Cia is given by Cia = h\fifj^ where h\ is a boson operator and fia is a fermion 
operator. In the slave fermion representation Cjo- = Instead of the difficult 

to handle constraint of Eg. Q , one considers more convenient slave-particle con¬ 
straints 


where the fermion (boson) operator keeps track of the spin and the boson (fermion) 
operator keeps track of the charge in the case of the slave-boson (slave-fermion) 
representation. Such slave-particle approaches are usually studied in a functional 
integral representation of the partition function with the no double occupancy 
constraints enforced with the help of Lagrange multiplier. To solve the problem the 
mean field approximation is usually applied and the Lagrange multiplier is taken 
to be independent of the lattice site. It means, however, that the local no double 
occupancy constraint is replaced by a global one with uncontrollable consequences. 


2. The Ising version of the t J model 

Because of the difficulties in solving the full t-J model, often its simplified versions 
are studied. One of them is the t-Jz- This model can be considered as a limiting 
(J_L = 0) case of the t-J model © which has an Ising rather than a Heisenberg 
spin interaction: 


Ht- 


t-J, 


t + Jz Sj ^nihj 

(bV (b> ^ 


(5) 


The original t-J model possesses the continuoues global SU(2) spin symmetry. In 
the Hamiltonian (|^ the interaction term SfSj has a lower discrete Z 2 symmetry. 
The rest of the terms, however, still possess the original SU(2) symmetry. As a 
result, the total symmetry of the t-Jz Hamilonian is dependend on the value of the 
Jz coupling. For = 0 the symmetry is SU(2) like in the full t-J Hamiltonian, 
whereas for = 0 it is only Z 2 . 
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Figure 1. Mapping between the physical states from the Hilbert space with no double occupied sites and 
the lattice spin and dopon states. 


In contradistiction to the t-J^ Hamilonian, all terms of the full Ising-t-J Hamil¬ 
tonian possess only the discrete Z 2 symmetry, independetly of the values of the 
model parameters. Unfortunately, since the operators cia transform themselves in 
the fundamental representation of SU(2), there is no obvious way to derive it from 
the t-J model ([^. One of the possibility is to use the enlarged spin-dopon repre¬ 
sentation of the operators Cia [TO! • 

In the framework of this approach fermion operators dia are assigned to doped 
carriers (like holes) instead to the lattice electrons. The vectors that span the 
enlarged on-site Hulbert space have the form of |(Ta), where a labels the 

2D lattice spin Qi Hilbert space and a = 0, t) i; labells the 4D onsite dopon 
Hilbert space. The physical subspace is spanned by the spin-up | fl 0)j, spin-down 
I 0)i, and spinless vacancy (| flDi — | /\/2 states [H]. The constraint 

QiMi + -nf = 0 , ( 6 ) 

has to be applied to remove the remaining unphysical states [12]. In the above 
Mr = E.,.' dl^Tfja'dia' is the dopon spin operator. This way the physical spin 
operator can be expressed as 


Si = Qi + M,. (7) 

Taking into accont that ((5“)^ = Eq. Q can be written as 

QfMt + nfY^ iQf)^ = 0. ( 8 ) 

a=x^y,z a=x,y,z 

In the full Ising t~J model the transverse spin components should vanish identically. 
According to Eq. Q, this requires Qf = Mf = 0. Then, the Ising t-J model can 
be derived by projecting the dopon operators onto the Hilbert space determined 
by the local constraint 

QfMf + = 0. (9) 

that is the Ising counterpart of Eq. (§. The projected physical electron operators 
Cia can be then expressed in terms of the lattice spin and dopon operators: 

= Q - 4^, (10a) 

= 'Pt4Tt = (2 + Ql) 4i 


(10b) 
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where the operator that projects onto the physical subspace is given by Vf^ = 
1 — (2Q?Mf + \nf). Then, it can be easily shown that the Ising version of Eq. 
is fulfilled: 


Sf = - 4c4) = Q! + Ml (11) 

and the transverse components of the physical spin operators vanish identically. 

s* = (s-y= 4^4 ^ O' (12) 

The adventage of the Ising representation of the t-J model is particularly visible 
close to half-filling, where the Hamiltonian ([^ is reduced to the following form 

IH: 


(bV (b> 


QfQ? - 7 


+ QfMj + 


(13) 


which has to be accompanied by Eq. (§. Since close to half-filling the hole con¬ 
centration 6 is small, in Eq. (13) we could drop the term describing the direct 
inter-dopon spin-spin interaction M?M?, which is proportional to 6'^. 

Since [Qf, = 0 the spin degrees of freedom in Eq. (13) can be described 

by classical variables. 

The constraint ([^ can be enforced with the help of a Lagrange multiplier. Since 
for each lattice site i QfM? + \nf > 0, the global Lagrange multiplier 




+ -Ai 
4 


(14) 


ensures that the constraint ([^ is fulhlled locally and the occupancy of an unphysical 
state at arbitrary site would lead to an increase of the total energy by A —)• -|-cx). 
As a result, all unphysical states are eliminated, so that the constraint ([^ fulfilled 
rigorously. 

The Hamiltonian © accompanied by the constraint ( |14[ ) represents a system 
described by classical (Q?) as well as quantum (d,) degrees of freedom. However, 
as pointed out above, the direct interaction between the quantum particles can 
be neglected close to half-filling and only the interaction between quantum and 
classical particles is present in Eq. (13). In this aspect, the Ising t~J model is 


similar to the Falicov-Kimball model and efficent hybrid methods that have been 
developed for latter model can be applied. 


3. Numerical approach 


The numerical technique we use to solve the effective model is based on a method 
that combines Monte Carlo simulations with exact diagonalization of one-particle 
Hamiltonians. This technique was proven to work effectively for the Falicov- 
Kimball model [IMS]. The details of the application of this method to the Ising 
t~J model are described in Ref. |10j . here we will sketch it for the sake of com¬ 
pleteness. 

The Hamiltonian of the Ising t-J model given by Eq. (13) can be divided into a 
one-particle part describing itinerant quantum particles with atomic levels varying 
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from site to site and a part describing Ising-type interactions between the classical 
variables Qf. The values of the atomic levels is determined by the distribution of 
the variables Qf. Together with the Lagrange multiplier term the Hamiltonian can 
be written as 


FaWdldj. + jY^ QIQ^^ + const. 




(U> 


(15) 


The hopping matrix Tija{^) is given by 




- + s{a)Qi 


{j)i 


(16) 


Since close to half filling details of the dispersion relation are important in strongly 
correlated systems [THHH] , in the above equation we used tij instead of the nearest- 
neighbour hopping t. By choosing proper values of tij = t{ri—rj) one can reproduce 
the dispersion relation of, e.g., high-Tc superconductors. In numerical calculations 
we restrict the hopping range to third nearest neighbours, i.e., only t, t' and t" 
are nonzero. s{a) is equal to 1 for a =f|' and -1 for a =1(; (j)j indicates that in the 
summation j runs over all nearest neighbours of site i. Note, that the quantum and 
classical degrees of freedom are coupled by the exchange constant J and by the 
Lagrange multiplier A. Numerical simulations indicate that both these couplings 
may be important, e.g, A is crucial in destroying an antiferromagnetic order when 
the concentration of holes increases, whereas J plays important role in formation 
of spin polarons. 


For a given distribution {Qf} of the classical variables the hopping matrix (16) 
can be numerically diagonalized and the Hamiltonian (15) can be rewritten as 


i/;i7(A) = ^({Qf},A)<X 


+ ^E«'3j 

(b> 


(17) 


where the constant term was neglected. This form of the Hamiltonian allows to 
carry out the classical Monte Carlo simulations based on a modified Metropolis 
algorithm [iSEn]. In the first step we choose an initial configuration {Qf}- It is 
defined by the distribution the lattice spins with three possibilities at each site: 
spin up, spin down, empty. The number of empty sites is given by the doping level. 
Depending on the physical problem some additional constraints may be imposed 
on the initial state. For example, we may require equal numbers of spin-up and 
spin-down sites to run a simulation in a subspace of the total magnetization equal 
to zero TiQi — Ti^^ — 0- Next, the Hamiltonian (15) is diagonalized and the 


free energy of the dopons in the initial state is calculated. Then, we attempt to 
change the configuration {Qf} —)• {Q'C- The changes can be twofold: one can 
modify the direction of one or two lattice spins or the distribution of the empty 
sites can be altered. The decision what kind of attempt is made is random. In the 
case of spin modifications if we work in a subspace of zero total magnetization, we 
randomly choose two lattice sites with opposite spin directions and exchange the 
spins. Otherwise we simply flip a randomly chosen spin. After the modification of 
the state is made, the Hamiltonian (15) is again diagonalized, what gives the en¬ 
ergy spectrum and the eigenstates of dopons. A new value of the dopon free energy 
is calculated and the configuration {Qf} is accepted or rejected according to the 
Metropolis criterion. This criterion is modified with respect to the original that is 
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used in simulations of classical systems: the internal energy in statistical weights 
is replaced by the free energy of the quantum subsystem (dopons). A detailed 
description of this approach can be found in Ref. m- The Monte Carlo simula¬ 
tion gives all the characteristics of both the classical (lattice spins) and quantum 
(dopons) subsystems, i.e., all correlation functions, specific heat, magnetization, 
spectral functions, etc. can be determined as a function of temperature, doping 
level, interaction strength, dispersion relation, etc. Moreover, since we work in the 
real space we can study inhomogeneous systems, e.g., with polarons. Since in this 
approach only a one-particle Hamiltonian has to be diagonalized there is no limit 
to the size of the system from the available computer memory. The only limit comes 
from the CPU time, because in each Monte Carlo step the matrix given by Eq. 


(16) is diagonalized, what significantly slows down the simulation in comparison 


to simulations of classical systems. Nevertheless, we are able to run simulations for 
50x50 lattices, what much beyond the capabilities of the fully quantum mechani¬ 
cal methods like the exact diagonalization based on the Lanczos algorithm or the 
Quantum Monte Carlo. 

The simulations have been carried out in the canonical ensemble, which allows 
for accurate control of the concentration of holes. The unphysical states have been 


removed by the term (14) with A of the order of a few hundreds. This way A is by 


far the largest energy scale in the system, which guaranties the single occupancy 
of each lattice site. 

One of the areas where the Ising t-J model can be applied is the problem of 
the rapid suppression of the antiferromagnetic order with increasing doping level 
in the high-Tc superconductors. In order to study the antiferromagnetic order we 
have to be able to calculate the spin-spin correlation function. In the Ising t-J 
model it can be defined as 


^ E E + M/)) j(r - |i^, - R,\), (18) 

i n 


where K = (vr, vr) and 


(I(x) = 


1 if \x\ < 0.5a, 
0 otherwise. 


with a being the lattice constant. (...) in Eq. (18) means an average over the 


spin configurations generated in the Monte Carlo run. This quantity will allow to 
describe the character of the antiferromagnetic correlations. For a long-range order 
it will has a finite value for arbitrary distance r, for a quasi-long-range order it will 
decay algebraically, and for a short-range order it will decay exponentially. Another 
quantity which is easy to calculate is the static spin-structure factor, given by 




(Iff) 


What is more interesting, this modified classical Monte Carlo approach can give also 
dynamic properties of the dopons, which are fully quantum mechanical particles. 
Namely, one can calculate the dopon’s spectral function 


A(k,(x) 


-Im G (k,(x + fO'*') , 

vr ^ ' 


( 20 ) 
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where 


G {k, z) = Yl e^HR.-R.) ( 

ija 


- - s{a)Q 


- - s{a)Q 


( 21 ) 


Here, similarly to Eqs. (18) and (19), (...) indicates averaging over spin configura¬ 
tions generated in the Monte Carlo runs and 


Q(7 (-^ 2 ? 


kl 



( 22 ) 


is the real-space Green function for a given spin conhguration {Qf} 
given by Eq. (16). Note that all the quantities given by Eqs. (18), 
are defined for physical electrons. 


Tklai^') is 
and dlol) 


4. Antiferromagnetism in the Ising t—J model 


The evolution of the antiferromagnetic Mott insulating state into a superconducting 
state is one of the most intriguing problems in the physics of the high-Tc super¬ 
conductors. In particular, it is difficult to explain how the antiferromagnetic order 
is destroyed very quickly when charge carriers are doped into a parent cuprate ma¬ 
terial. In most thermodynamic measurements for hole doped cuprates, long range 
antiferromagnetism does not coexist with superconductivity and disappears com¬ 
pletely around doping density 6 ~ 5%. Most of analytical and numerical studies 
of the t-J model show that while upon doping the antiferromagnetic order is sup¬ 
pressed, it survives to much larger hole density that observed in experiments. The 
discrepancies may imply that the t-J model is insufficient to describe the physics 
of the hight-Tc superconductors. But they may also imply that the methods used 
to study this model close to half filling are not reliable enough to give the correct 
value of the critical hole density. It was already mentioned in the Introduction that 
most of both analytical and numerical methods have difficulties in dealing with 
the t-J model close to half filling, where the constraint of no double occupancy 
is particularly important. This is the regime where we believe the validity of the 
Ising version of the t-J model is most justified. 

Most of the results for the t-J model close to half filling are restricted to one or 
two holes in an antiferromagnetic background |211430) . These methods do take into 
account the strong electron correlations, however, they do not allow to change the 
hole concentration and study the evolution of the antiferromagnetic order. On the 
other hand, variations of mean-field-type approaches 1314440] allow to control the 
hole density, but their validity is questionable in the underdoped regime, where 
electronic correlations are crucial due to the proximity to the Mott state. 

The proposed numerical approach to the Ising t-J model takes advantages from 
both these groups of methods: on the one hand the concentration of holes can 
be changed almost continuously from zero to an arbitrary density, on the other 
hand the no-double-occupancy constraint is fulfilled rigorously not on average, 
like in the mean-field approaches, but at every lattice site. This was possible at 


the expense of neglecting the transverse spin-flip term (12). However, it was shown 
in Refs. [TO] and mi that the energy of one and two holes calculated for the full 
t-J model |424444j and for its Ising version are close. 

Eigure shows the spin-spin correlation function g{r) defined by Eq. 18 for 
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Figure 2. Distance dependence of the spin-spin correlation function g{r) for different hole contrentations 
5 = 0.02 (a), 0.04 (b), 0.06 (c) and 0.08 (d). A logarithmic scale is used on the vertical axis. The solid lines 
show fits to the Monte Carlo results. The following parameters have been assumed: J = 0.2t, kT = O.lt, 
t' = —0.27t and t" = 0.2t. 


different hole concentrations. In order to describe the decay of the correlations we 
use a logarithmic scale on the vertical axis. One can see in this figure that the 
character of this correlation changes very rapidly with the increase of the number 
of holes. For a very small concentration 6 = 0.02 (Fig. [^) the correlation drops 
at a very short distance but then it is almost constant for larger distances, what 
indicates the presence of the long range antiferromagnetic ordei[^ For a slightly 
higher doping 6 = 0.04 (Fig. [^) we can observe an exponential decay at small 
distance, but then it slows down at larger distance changing into an algebraic 
decay. It suggests the presence of the quasi-long-range order. Finally, when we 
further increase the hole concentration to <5 = 0.06 and 6 = 0.08 (Figs. and 
Hi, respectively), we observe an exponential decay at all distances, what means 
that the long range antiferromagnetic has been destroyed. It may suggest that 
the critical hole concentration in the Ising t-J model is below 6%, what is in a 
perfect agreement with experiments. This result, however, has been obtained on a 
relatively small cluster and should be confirmed by the finite-size scaling. 

There is still, however, the question about the nature of the suppression of the 
long range antiferromagnetic order. Writing explicitly the A-dependent term in 
Eqs. @ and @ 




l + Q 


+ (^2 



(23) 


^The presented results were calculated for a 20x20 system with periodic boundary conditions and therefore 
we cannot say anything about the behaviour of g{r) at a very large distance. 
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5 = 0.02 S = 0.04 

a) b) 



S = 0.12 




Figure 3. Spectral functions for different hole concentration. 


one can see that in a perfect Neel state single holes which hop only to nearest 
neighbours are fully localized. The holes can gain kinetic energy by destruction of 
the antiferromagnetic order and forming a ferromagnetic region (spin polaron), but 
this mechanism is effective only for a very small value of the exchange J [411 05] . 
Nonzero values of t' and t" allow holes to propagate and gain some energy by intra¬ 
sublattice hoppings. The same situation persists for a small but finite concentration 
of holes. The hole spectral function calculated according to Eq. (20) for the hole 
concentration <5 = 0.02 is presented in Fig. [^. It is exactly the spectral function 
for free electrons with the dispersion relation given by hoppings only to the second 
and third neighbours (i.e., within the same sublattice) with the hopping integrals t' 
and t", respectively. When the hole concentration increases the potential gain from 
their mobility enhancement would increase as well. And at some point it starts to 
be energetically favourable to locally destroy the antiferromagnetic order and allow 
the holes to hop also to nearest neighbour sites. The evolution of this process can 
be observed in Figs.[^-f, where the contribution from electrons with the dispersion 
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relation characteristic for nearest neighbour hopping becomes more and more pro¬ 
nounced. This results suggest that the main mechanism that is responsible for the 
suppression of the long-range antiferromagnetic order is the competition between 
the hole mobility and tendency towards minimization of the spin-spin exchange 
energy. This strong competition results directly from the no double occupancy con¬ 
straint. This constraint, in connection with the lack of the spin-flip term in the 
Ising t-J model prevents holes from hopping to nearest neighbour sites. If the gain 
is comparable to the antiferromagnetic ordering energy the total energy can be 
lowered by suppressing the order. 


5. Nagaoka polar on 

Since the proposed numerical approach to the Ising t-J model does not require 
translational invariance it is well suited to study inhomogeneous systems. This 
feature allows to study the small-J and small hole concentration limits of the 
t-J model. In the small-J limit the dynamics of spins is much slower then the 
dynamics of charge carriers what may justify the approximation that leads to the 
Ising version of the t-J model, e.g., neglecting of the transverse spin components 
[Eq. @]. 

According to the Nagaoka theorem [451 00] the infinite-!/ Hubbard model with 
only a single hole has a fully spin-polarized ferromagnetic ground state. This regime 
corresponds to the t~J model in the J —?■ 0 limit. The Nagaoka theorem is excep¬ 
tional in the sense that it is one of very few rigorous results for strongly correlated 
systems. Unfortunately, its validity is limited only to a one particular case. The 
situation for large but finite U (what is equivalent to J > 0) and/or small but finite 
hole concentration is much less clear. For finite J there is the antiferromagnetic 
exchange energy that competes with the kinetic energy of holes in a ferromagnetic 
spin background. Therefore a single hole in a system with finite J may lead to 
formation of a ferromagnetic ’’bubble” that allows the hole to gain the kinetic 
energy. The rest of the system would have anitferromagnetically ordered spins to 
minimize the exchange energy. The size of the ferromagnetic region would be de¬ 
termined by the competition between the exchange and kinetic energies in a way 
that guarantees a global minimum of the total |42) . 

The situation becomes more complicated with the increase of the number of holes. 
Already for two holes results are ambiguous. In the small-J regime the question 
is whether the two holes will form a single bipolaron or two sparate polarons. 
Or more generally, whether the holes form a bound state. The problem occurs 
because the characteristic length scale in the small-J regime, connected with the 
size of the ferromagnetic region, is large, beyond the limits of applicability of most 
of the fully quantum-mechanical approaches. One of the few methods which are 
capable to calculate properties of two holes in an antiferrormagnetic background is 
an accurate exact diagonalization method, defined over a limited functional space 
(EDLFS) recently proposed by Bonca et al. in Refs. [l3llll]. The construction of 
the limited space starts from a Neel state with two holes located on neighboring 
lattice sites. Next, the kinetic and the spin-flip parts of the t-J Hamiltonian are 
applied to generate the basis vectors. The ground state is then calculated within the 
generated functional space by means of the Lanczos method. This approach allows 
to study much larger systems than the standard exact diagonalization methods. 
Nevertheless, the maximum distance between the holes is limited by the size of 
the generated functional space. The problem is that the size of the ferromagnetic 
region which may contain the two holes diverges with decreasing J and at some 
point even this method does not allow to study the bipolaron problem. The value 
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Figure 4. Snapshot of the lattice spin configurations and false color plots of corresponding two lowest hole 
wave functions for J = 0.01. In the leftmost picture filled (blue) circles represent lattice spins pointing up 
and the empty (white) ones lattice spins pointing down. The white rectangular represent a ferromagnetic 
spin polaron. 


of J below which the EDLFS method cannot give reliable results is about 0.04. On 
the other hand, the approximations that lead to the Ising t-J model can be applied 
for arbitrarily small J with its accuracy increasing with decreasing J. Therefore, 
we used a comparison of the results for the Ising t-J model and the results of the 
EDLFS method applied to the full t-J model to examine the validity of neglecting 
the transverse spin components and to get insight into the physical meaning of 
this approximation. One important parameter that can be compared is the average 
distance between two holes. In order to calculate its value in the Ising t-J model we 
run Monte Carlo simulations for a system with two holes. The difference between 
the present simulations and those carried out for the study of the destruction of 
the antiferomagnetic order with the increase of the hole concentration is that now 
we do not keep zero total magnetization. In the previous simulations each Monte 
Carlo attempt consisted of two spin flips of two opposite lattice spins. Here the have 
two types of attempts: a transfer of a randomly chosen lattice spin from one site to 
another or a single spin flip. The latter kind of attempts does not conserve the total 
lattice magnetization. Figure shows a typical low-temperature configuration of 
the lattice spins and two lowest corresponding hole wave functions Ti and T 2 . The 
average distance between the holes is calculated as 

«=(ee \ri-r2\V{ri,r2)\ , (24) 

\ ri r2 / 


where ri, r 2 run over all lattice sites, 


V{ri,r2) 


'^2{ri) '^2{r2) 


(25) 


and (...) denotes an average over the lattice spin configurations generated in a 
Monte Carlo run. It turned out that the distance between two holes in the Ising 
t-J model has the same dependence on J as in the full t-J model, but the numeric 
prefactor is almost 30% smaller. Namely, the distance D{J) in the Ising t-J model 
is given by 1.4 and 1.97 in the full t-J model. The letter function has 

been obtained from a finite size scaling of the results of the EDLFS method m- 
The Monte Carlo simulations were carried out for J < 0.04 and the EDLFS method 
was used for J > 0.04. For such a small value of J the Monte Carlo results indicated 
that the energy of one bipolaron is smaller than that of two polarons, what suggests 
binding of the holes. The difference between the hole-hole distance in the full t-J 
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J = 0.01 J = 0.03 J = 0.05 J = 0.07 J = 0.09 J = 0.11 J = 0.13 J = 0.15 



Figure 5. Snapshots of the lattice spin configurations for different numbers of holes and for different 
exchange interaction J. The meaning of the small filled and empty circles is the same as in Fig. The 
main polaron is formed by lattice spins pointing down, whereas the dark areas in configuration for J > 0.07 
and 8 and 10 holes represent separate spin—up ferromagnetic polarons. 


model and in its Ising version can be explained by the approximations used in the 
Ising t-J model: On the one hand, when the transverse components of the spin 
operators are neglected the boundary between the ferromagnetic ’’bubble” and its 
antiferromagnetic surroundings is impenetrable. On the other hand, the spin flip 
term in the full t-J model allows a hole to enter the antiferromagnetic region. 
The movement of a hole through this region is accompanied by formation of a 
string of defects in the antiferromagnetic order what strongly limits the range this 
penetration. The behaviour of the holes in these two models can be explained with 
the help of an analogy to particles in quantum wells: the case of the Ising t-J model 
it would be an infinite rectangular quantum well, whereas in the full t-J model the 
walls of the well would be inclined outward, what would lead to a slightly broader 
wave function. 


5.1. Finite density of holes 

For two holes in the small~J limit it is energetically favourable to form a single 
ferromagnetic spin where the holes can move freely. Then, the question is whether 
this scenario will hold for higher number of holes. Figurej^show snapshots of Monte 
Carlo simulations for up to 10 holes in a 20x20 system for J from 0.01 to 0.15. 
One can see there that for J > 0.07 and 8 and 10 holes multiple polarons are 
formed. This is, however, the region where the spin-flip processes may play a more 
important role and the validity of the approach may be questionable. Therefore, we 
restrict ourselves to J < 0.05, similarly to the case of two holes. The Monte Carlo 
studies in this regime show that the energy E as a function of the number of holes 
can be well fitted by aN+by/N, where b is positive. It means that the function E{N) 
is concave and for all the studied hole concentrations (N < 10) it is energetically 
favourable to phase separate the system in a hole-rich ferromagnetic region and 
an antiferromagnetic region without holes. As can be seen in Fig.[^ the size of the 
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Figure 6. Fraction of the total number of the lattice sites occupied by the ferromagnetic polaron as a 
function of J for different holes numbers. 

ferromagnetic ’’bubble” decreses with increasing J what can be explained by the 
increasing cost of broken antiferromagnetic bonds within the spin polaron. Fig. 
shows the size of the ferromagnetic polaron as a function of J. The points from 
Monte Carlo simulations are there fitted by a function Ns{J)/N = aF, where Ng is 
the number number of the lattice sites occupied by the ferromagnetic polaron, N is 
the total number of the lattice sites and a and b are fitting parameters. For a fixed 
J the size of the polaron increases with increasing hole concentration and at some 
point it includes all the lattice sites. This situation resembles the Nagaoka state, 
but for a finite number of holes. With decreasing J the critical hole concentration 
decreases and in the J —)• 0 limit {U —>■ oo) this state is converted into the standard 
Nagaoka state with vanishing hole concentration. The dependence of the critical 
hole concentration as a function of J can be well fitted by 5t = 0.44 what is 
very close to where the latter form can be easily derived by comparing the 

kinetic energy of a few holes in an otherwise empty band and the exchange energy of 
the lattice spins. If the hole concentration exceeds 5t all the lattice spins are fully 
polarized. For a hole concentration lower than 6t the system is phase separated 
into a hole-rich ferromagnetic part and a hole-depleted antiferromagnetic part. In 
this regime the size of the hole-rich ferromagnetic polaron (the number of lattice 
sites) linearly depends on the hole concentration. The above results hold true in 
the limit of small J. As can be seen in Fig. for J > 0.07 holes are confined to 
several separate polarons. In this regime, however, the validity of neglecting the 
transverse spin components is questionable and the multi-polaron picture may not 
be relevant to the isotropic t~J model. 


6. Summary 

We have presented a representation of the of the t~J model where the system is 
described in terms of fermions interacting with static localized spins. Although it is 
a slave-particle approach, in contrast with many similar approaches, the local no- 
double-occupancy constraint is rigorously taken into account. Within the proposed 
approach we have shown that the long range antiferromagnetic order disappears 
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already at the doping of the order of a few percent, what is in an agreement with 
the experimental data for high-Tc superconductors. Additionally, we have demon¬ 
strated that the no-double-occupancy constraint is responsible for the destruction 
of the order. Since it is difficult to take this constraint into account in most of the 
analytical and numerical approaches to the t-J model, it explain why theoretical 
estimations of the critical dopant concentration usually give a significantly larger 
value. 

The proposed approach to the t~J model does not require translational invariance 
of the system, what allow to use it to study inhomogeneous systems. Exploiting 
this feature we have studied also formation of the Nagaoka polaron in the small-J 
limit. The main difficulty in analysing the t-J model in this limit is that a large 
size of the lattice is required to correctly describe the dynamics of holes. In the 
proposed approach, however, the lattice spins are treated as classical variables, 
what allows to study systems much larger than in fully quantum approaches like 
the Lanczos method, quantum Monte Carlo, DMRG, EDLFS, etc. Therefore, we 
were able to show that it is energetically favourable for the system to segregate 
into the ferromagnetic hole-rich phase and hole-depleted antiferromagnetic phase. 
The size (surface) of the ferromagnetic bubble depends linearly on the number of 
holes, while its dependence on J is given by the square-root function. 
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